/**
 * @file PRE79Simulation.h
 *
 * Simulation classes
 *
 */

#ifndef _PRE79SIMULATION_H
#define _PRE79SIMULATION_H

#include "PRE79StandardHamiltonian.h"
#include "Simulation.h"
#include "serializer.h"
#include "ILoggable.h"
#include "SimulationDBFind.h"
#include "AutoCorrelationTimeCalculator.h"
#ifndef OMP_DEBUG
#include <omp.h>
#else
#include "omp_debug.h"
#endif

/**
 * A redundant class for a single production when a parallel production is performed.
 */
class PRE79Production: public ILoggable
{
    unique_ptr<Lattice>                     lattice;            ///<pointer to the lattice object
    unique_ptr<StandardL2Hamiltonian>    H;                  ///<pointer to the Hamiltonian object
    unique_ptr<LatticeSimulation>           simulation;         ///<pointer to the LatticeSimulation object
    unique_ptr<PRE79StandardProperties>     prop;               ///<pointer to the properties calculator and accumulator
    unique_ptr<Metropolis>                  metro;              ///<pointer to the Metropolis method
    const Settings *   settings;     ///<local copy of settings
    long nprod;                 ///<total number of production cycles
    long ncycles;               ///<number of cycles considered for measurement
public:
    /**
     * Constructor
     *
     * @param set       reference to global Settings object
     * @param _nprod    number of production cycles for the entire production
     * @param start     initial state for production
     */
    PRE79Production(const Settings & set, long _nprod, const Lattice & start):
        nprod(_nprod),
        settings(&set)
    {
        ncycles = nprod / settings->simulation.measure_frequency;

        ///create the needed objects
        lattice = make_unique<Lattice>(start);
        H = make_unique<StandardL2Hamiltonian>(settings->hamiltonian.temperature, settings->hamiltonian.lambda, settings->hamiltonian.tau, settings->hamiltonian.h, settings->hamiltonian.kappa);
        metro = make_unique<Metropolis>(*settings, H.get());
        prop = make_unique<PRE79StandardProperties>(lattice.get(), ncycles);
        simulation = make_unique<LatticeSimulation>(H.get(), lattice.get(), metro.get(), nprod, 0);
        SetStream(&std::cout);
    }
    /// disable copying
    PRE79Production(const PRE79Production &) = delete;

    /**
     * Run a single production in the set.
     */
    void Run()
    {

        #pragma omp critical
        Log() << "Production with freq " << settings->simulation.measure_frequency << std::endl ;

        ///save the time for calculation of remaining time
        pt::ptime start_t = pt::second_clock::local_time();
        ///the remaining time will be reported 5 times through the production
        int remaining_interval = simulation->GetNCycles() / 5;
        if(settings->simulation.calculate_time)
            Log() << "Remaining time will be reported every " << remaining_interval << " cycles\n";

        ///an instance of the autocorrelation time calculator
        AutoCorrelationTimeCalculator ac(
            lattice.get(),
            settings->simulation.autocorrelation_frequency,
            settings->simulation.autocorrelation_length);

        ///k will count the cycles
        long k = 0;
        ///this loop will terminate if the underlying #LatticeSimulation Iterate()
        ///funcion returns false, which will happen when all the cycles are finished
        ///the Iterate() function does the actual lattice sweeps
        while(simulation->Iterate())
        {
            ///update the autocorrelation time calculator
            ac.Update();

            ///do measurements here
            if(k % settings->simulation.measure_frequency == 0)
            {
                ///measurements are done via the Update() function of the properties object
                prop->Update(k, H.get(), &ac);
            }

            ///here the progress as percentage is calculated and reported
            if(k % 1000 == 0 && settings->output.report_progress)
            {
                Log() << "E = " << prop->EnergyEvolution()[k / 1001] << std::endl;
                Log() << "Progress: " << (double(k) / double(simulation->GetNCycles())) * 100.0 << "%\n";
            }


            ///runtime adjustment of the random-walk radius
            // Swendson (2011) the interval shouldn't be less than (total number of MC cycles)^(1/2)
            if(k % settings->simulation.radius_adjustment_frequency == 0)
                metro->AdjustRadius(lattice.get());

            ///measure the acceptance rate if needed
            if(settings->simulation.measure_acceptance && k % settings->simulation.measure_acceptance_frequency == 0)
            {
                //Log() << "Acceptance rate: " << metro->MeasureAccepted(lattice)*100.0 << "%\n";
                Log() << "Acceptance rate: " << simulation->GetAcceptance() * 100.0 << "%\n";
                Log() << "Mean acceptance rate: " << simulation->GetMeanAcceptance() * 100.0 << "%\n";
            }

            ///calculate the remaining time of the present production
            if((k + 1) % remaining_interval == 0 && settings->simulation.calculate_time)
            {
                pt::time_duration run1k = pt::second_clock::local_time() - start_t;
                int total1kruns = (simulation->GetNCycles() - k - 1) / remaining_interval;
                std::cout << "Thread: " << omp_get_thread_num() << "/" << omp_get_num_threads() <<  ": " << pt::to_simple_string(run1k * total1kruns) << " remaining\n";
                start_t = pt::second_clock::local_time();
            }
            k++;
        }
    }

    ///Access the pointer to the lattice object
    const Lattice * GetLattice() const
    {
        return lattice.get();
    }
    ///Access the pointer to the properties object
    const PRE79StandardProperties * GetProperties() const
    {
        return prop.get();
    }
    ///Access the pointer to the LatticeSimulation object (read only)
    const LatticeSimulation * GetSimulation() const
    {
        return simulation.get();
    }
    /**
     * Redefinition of the Log stream, which adds thread information to the output.
     * The implementation at #ILoggable doesn't have this.
     */
    virtual loggable_proxy<ILoggable> & Log()
    {
        ILoggable::Log() << "Thread: " << omp_get_thread_num() << "/" << omp_get_num_threads() << ": ";
        return ILoggable::Log();
    }

};

/**
 * This class takes care of the entire process of preparing the simulation, thermalizing the initial states
 * and the actual production phase. Several mechanisms are implemented, such as
 * searching for thermalized states.
 */
class PRE79Simulation: public ILoggable
{
    unique_ptr<Lattice>                     lattice;            ///<pointer to the lattice instance
    unique_ptr<StandardL2Hamiltonian>    H;                  ///<pointer to the hamiltonian object
    unique_ptr<LatticeSimulation>           thermalization;     ///<underlying thermalization phase
    unique_ptr<LatticeSimulation>           simulation;         ///<underlying production phase
    unique_ptr<PRE79StandardProperties>     prop;               ///<pointer to production properties calculator
    unique_ptr<Metropolis>                  metro;              ///<pointer to the Metropolis method
    std::vector<unique_ptr<PRE79Production>>    productions;   ///<a vector of productions if parallel production is enabled
    unique_ptr<PRE79StandardProperties>     thermalprops;       ///<pointer to thermalization properties calculator
    const Settings  *   settings;     ///<local settings reference
    SimulationDB        database;       ///<simulation database interface
    bool    restored;      ///<should we search for a thermalized state to reuse as an initial state?

    /**
     * Initialization
     */
    void Init()
    {
        Log() << "Creating Hamiltonian\n";
        H = make_unique<StandardL2Hamiltonian>(settings->hamiltonian.temperature, settings->hamiltonian.lambda, settings->hamiltonian.tau, settings->hamiltonian.h, settings->hamiltonian.kappa);
        Log() << "Creating Metropolis\n";
        metro = make_unique<Metropolis>(*settings, H.get());

        if(!restored)
        {

            ///--- search for a possible somewhat-thermalized state
            if(settings->simulation.find_thermalized)
            {
                Log() << "Searching database for thermalized state\n";
                bool found = false;

                int cycle = 0; ///< will store the cycle at which the state was left, unused

                ///--- allow for a temperature tolerance in searching for states (subsequent thermalization is still necessary)
                Lattice found_state;
                if(settings->simulation.find_thermalized_temperature_tolerance != 0.0)
                {
                    Log() << "Temperature tolerance enabled\n";
                    found_state = FindLastStateTemperatureTolerant(*settings, found, cycle, settings->simulation.find_thermalized_temperature_tolerance);
                }
                else
                {
                    found_state = FindLastState(*settings, found, cycle);
                }
                ///---
                ///--- allow for a field tolerance in searching for states (subsequent thermalization is still necessary)
                if(settings->simulation.find_thermalized_h_tolerance != 0.0)
                {
                    Log() << "Field tolerance enabled\n";
                    found_state = FindLastStateFieldTolerant(*settings, found, cycle, settings->simulation.find_thermalized_h_tolerance);
                }
                else
                {
                    found_state = FindLastState(*settings, found, cycle);
                }
                ///---

                if(found)
                {
                    lattice = make_unique<Lattice>(found_state);

                    Log() << "Creating Supplementary Thermalization for recovered state\n";
                    thermalization = make_unique<LatticeSimulation>(H.get(), lattice.get(), metro.get(), settings->simulation.supplementary_thermalization_cycles);
                }
                else
                {
                    Log() << "No thermalized state found, creating Thermalization\n";
                    thermalization = make_unique<LatticeSimulation>(H.get(), lattice.get(), metro.get(), settings->simulation.thermalization_cycles);
                }
            }
            //---
            else
            {
                Log() << "No thermalized state found, creating Thermalization\n";
                thermalization = make_unique<LatticeSimulation>(H.get(), lattice.get(), metro.get(), settings->simulation.thermalization_cycles);
            }
        }
        else
        {
            Log() << "Previous final state passed on as initial, creating Supplementary Thermalization\n";
            thermalization = make_unique<LatticeSimulation>(H.get(), lattice.get(), metro.get(), settings->simulation.supplementary_thermalization_cycles);
        }
        ///--- UNFINISHED - resume simulation from a saved state
        // TODO: the state is loaded, but the runtime properties are not present and need to be calculated or loaded...
        int found_cycle = 0;
        if(settings->simulation.pick_up_aborted)
        {
            Log() << "Searching database for aborted simulation\n";

            bool found = false;
            Lattice found_state = FindLastState(*settings, found, found_cycle);
            if(found)
            {

                lattice = make_unique<Lattice>(found_state);
                Log() << "Recovering from aborted simulation at " << found_cycle << " cycles\n";
            }
            else
            {
                Log() << "No aborted simulation found\n";
            }
        }
        //---
        int cycle_advantage = settings->simulation.production_cycles - found_cycle;

        Log() << "Creating Production\n";
        simulation = make_unique<LatticeSimulation>(H.get(), lattice.get(), metro.get(), cycle_advantage, found_cycle);
        Log() << "Creating Properties\n";
        prop = make_unique<PRE79StandardProperties>(lattice.get(), settings->simulation.production_cycles / settings->simulation.measure_frequency);
        thermalprops = make_unique<PRE79StandardProperties>(lattice.get(), thermalization->GetNCycles() / 10);

    }

public:
    PRE79Simulation(const Settings & set):
        settings(&set),
        database(set),
        thermalprops(nullptr), thermalization(nullptr), lattice(nullptr)
    {
        if(set.initial.isotropic)
            lattice = make_unique<Lattice>(set, Lattice::Isotropic);
        else if(set.initial.biaxial)
        {
            if(set.initial.righthanded)
                lattice = make_unique<Lattice>(set, Lattice::BiaxialRighthanded);
            else
                lattice = make_unique<Lattice>(set, Lattice::Biaxial);
        }
        else if(set.initial.biaxial_alt)
        {
            if(set.initial.righthanded)
                lattice = make_unique<Lattice>(set, Lattice::BiaxialRighthandedAlt);
            else
                lattice = make_unique<Lattice>(set, Lattice::BiaxialAlt);
        }
        else if(set.initial.righthanded)
            lattice = make_unique<Lattice>(set, Lattice::IsotropicRighthanded);
        else
            lattice = make_unique<Lattice>(set, Lattice::Isotropic);
        restored = false;
        Init();
    }

    /// new simulation from a saved state
    PRE79Simulation(const Settings & set, const Lattice & saved):
        settings(&set), database(set)
    {
        restored = true;
        lattice = make_unique<Lattice>(saved);
        Init();
    }
    PRE79Simulation(const PRE79Simulation &) = delete;

    /// measure the duration of a 1000 cycles (for time estimation)
    pt::time_duration DurationOf1000Cycles() const
    {
        LatticeSimulation test(H.get(), lattice.get(), metro.get(), 1000);
        PRE79StandardProperties testprop(lattice.get(), 1000 / settings->simulation.measure_frequency);
        pt::ptime start_t = pt::second_clock::local_time();

        int k = 0;
        while(test.Iterate())
        {
            if(k % settings->simulation.measure_frequency == 0)
                testprop.Update(k, H.get());
            k++;
        }
        pt::time_duration run1k = pt::second_clock::local_time() - start_t;
        return run1k;
    }
    pt::time_duration ExpectedSimulationTime() const
    {
        int total1kruns = (thermalization->GetNCycles() + simulation->GetNCycles()) / 1000 ;
        return DurationOf1000Cycles() * total1kruns;
    }

    /// thermalization
    const Lattice & Thermalize()
    {
        AutoCorrelationTimeCalculator ac(lattice.get(), settings->simulation.autocorrelation_frequency, settings->simulation.autocorrelation_length);
        Log() << "Thermalization cycles: " << thermalization->GetNCycles() << std::endl;
        int tcycle = 0;
        while(thermalization->Iterate())
        {
            ac.Update();
            if(tcycle % 10 == 0)
                thermalprops->Update(tcycle, H.get(), &ac);
            if(tcycle % 1000 == 0 && settings->output.report_progress)
            {
                Log() << "E = " << thermalprops->EnergyEvolution()[tcycle / 1001] << std::endl;
                Log() << "Progress: " << (double(tcycle) / double(thermalization->GetNCycles())) * 100.0 << "%\n";
            }
            ///--- CAUTION - live adjustment of the random walk radius - could break detailed balance!!!
            if(tcycle % settings->simulation.radius_adjustment_frequency == 0)
                metro->AdjustRadius(lattice.get());
            ///---
            if(settings->simulation.measure_acceptance && tcycle % settings->simulation.measure_acceptance_frequency == 0)
            {
                Log() << "Acceptance rate: " << thermalization->GetAcceptance() * 100.0 << "%\n";
                Log() << "Mean acceptance rate: " << thermalization->GetMeanAcceptance() * 100.0 << "%\n";
            }
            tcycle++;
        }
        if(settings->output.save_thermalization_properties)
        {
            Log() << "Saving thermalization history\n";
            database.StoreThermalizationHistory(*settings, *thermalprops);
        }

        return *lattice;
    }

    ///
    /// One single-threaded thermalization followed by multi-threaded production
    /// (replicas). The replicas start from the same thermalized state, but
    /// use different seeds for the random number generator.
    /// \return The final state with minimum energy
    ///
    Lattice RunParallel()
    {
        ///--- thermalization
        Log() << "Adjusting radius\n";
        metro->AdjustRadius(lattice.get());
        if(settings->simulation.calculate_time)
        {
            Log() << "Calculating expected duration of simulation\n";
            pt::time_duration   expected = DurationOf1000Cycles() * (thermalization->GetNCycles() + settings->simulation.production_cycles / settings->openmp.number_of_threads) / 1000;
            Log() << "Expected time of simulation: " << pt::to_simple_string(expected) << std::endl;
        }
        if(settings->simulation.thermalization_cycles > 0 && settings->simulation.supplementary_thermalization_cycles > 0)
        {
            Log() << "Thermalization\n";
            Thermalize();
        }
        ///---

        ///--- set up partial simulations (on replicas)

        /// here there can be no doubt as for the number of currently running threads
        omp_set_dynamic(0);
        omp_set_num_threads(settings->openmp.number_of_threads);

        for(int i = 0; i < settings->openmp.number_of_threads; i++)
            productions.push_back(make_unique<PRE79Production>(*settings, settings->simulation.production_cycles / settings->openmp.number_of_threads, *lattice));

        Log() << "Starting " << settings->openmp.number_of_threads << " productions\n";

        /// GO!
        #pragma omp parallel for schedule(runtime) shared(rng2) private(random01)
        for(int i = 0; i < productions.size(); i++)
        {
            productions[i]->SetStream(this);
            productions[i]->Run();
        }

        ///
        /// Concatenate the resulting properties
        ///
        PRE79StandardProperties generalprop(lattice.get(),0);
        for(int i = 0; i < productions.size(); i++)
        {
            generalprop.Append(productions[i]->GetProperties());
        }
        generalprop.CalculateSpecificHeat();

        if(settings->output.save_properties_evolution)
        {
            Log() << "Saving Properties Evolution\n";
            database.StorePropertiesEvolution(*settings, generalprop);
        }
        if(settings->output.save_final_configuration)
        {
            Log() << "Saving Final Lattice\n";
            database.StoreFinalLattice(*settings, *productions[0]->GetLattice());
        }
        if(settings->output.save_final_properties)
        {
            Log() << "Saving Final Properties\n";
            PRE79MeanProperties pp(&generalprop, H.get());
            database.StoreFinalProperties(*settings, pp);
        }

        Log() << "Done\n";
        Log() << "Mean EPM: " << generalprop.TemporalMeanEnergyPerMolecule().Print() << std::endl;
        Log() << "Specific Heat: " << generalprop.SpecificHeat().Print() << std::endl;

        ///return the state with the lowest energy
        vect weights(0.0, productions.size());
        for(int i = 0; i < weights.size(); i++)
        {
            weights[i] = productions[i]->GetLattice()->GetMeanEPM();
        }
        int best = MinimumIndex(weights);

        Lattice ret = *productions[best]->GetLattice();
        ///--

        productions.clear();

        return ret;
    }


    ///
    /// Single-threaded simulation.
    /// \return Resulting state
    ///
    const Lattice & Run()
    {
        Log() << "Production cycles: " << simulation->GetNCycles() << std::endl;
        Log() << "Adjusting radius\n";
        metro->AdjustRadius(lattice.get());
        if(settings->simulation.calculate_time)
        {
            Log() << "Calculating expected duration of simulation\n";
            pt::time_duration   expected = ExpectedSimulationTime();
            Log() << "Expected time of simulation: " << pt::to_simple_string(expected) << std::endl;
        }
        if(settings->simulation.thermalization_cycles > 0 && settings->simulation.supplementary_thermalization_cycles > 0)
        {
            Log() << "Thermalization\n";
            Thermalize();
        }

        ///--- intermediate state saving frequency
        int intermediate_frequency = simulation->GetNCycles() / settings->output.intermediate_states;
        ///---

        #pragma omp critical
        Log() << "Production with freq " << settings->simulation.measure_frequency << std::endl ;
        long k = 0;

        pt::ptime start_t = pt::second_clock::local_time();
        int remaining_interval = simulation->GetNCycles() / 20 + 1;
        if(settings->simulation.calculate_time)
            Log() << "Remaining time will be reported every " << remaining_interval << " cycles\n";

        AutoCorrelationTimeCalculator ac(lattice.get(), settings->simulation.autocorrelation_frequency, settings->simulation.autocorrelation_length);

        while(simulation->Iterate())
        {
            ac.Update();
            ///--- update the statistical properties and store lattice state
            if(k % settings->simulation.measure_frequency == 0)
            {
                prop->Update(k, H.get(), &ac);
                if(settings->output.save_configuration_evolution)
                {
                    database.StoreLattice(*settings, *lattice, k);
                }
            }
            ///---

            ///--- progress report
            if(k % 1000 == 0 && settings->output.report_progress)
            {
                Log() << "E = " << prop->EnergyEvolution()[k / 1001] << std::endl;
                Log() << "Progress: " << (double(k) / double(simulation->GetNCycles())) * 100.0 << "%\n";
            }
            ///---

            ///--- intermediate snapshots of average statistical properties
            if(k % intermediate_frequency == 0)
            {
                if(settings->output.save_intermediate_states)
                {
                    database.StoreLattice(*settings, *lattice, k);
                    PRE79MeanProperties pp(prop.get(), H.get());
                    database.StoreProperties(*settings, pp, k);
                }

            }

            ///--- CAUTION: excessive adjustment of Monte Carlo step radius breaks detailed balance
            // Swendson (2011)
            if(k % settings->simulation.radius_adjustment_frequency == 0)
                metro->AdjustRadius(lattice.get());
            ///---

            ///--- measurement of acceptance rate
            if(settings->simulation.measure_acceptance && k % settings->simulation.measure_acceptance_frequency == 0)
            {
                Log() << "Acceptance rate: " << simulation->GetAcceptance() * 100.0 << "%\n";
                Log() << "Mean acceptance rate: " << simulation->GetMeanAcceptance() * 100.0 << "%\n";
            }
            ///---

            if((k + 1) % remaining_interval == 0 && settings->simulation.calculate_time)
            {
                pt::time_duration run1k = pt::second_clock::local_time() - start_t;
                int total1kruns = (simulation->GetNCycles() - k - 1) / remaining_interval;
                std::cout << "Thread: " << omp_get_thread_num() << "/" << omp_get_num_threads() <<  ": " << pt::to_simple_string(run1k * total1kruns) << " remaining\n";
                start_t = pt::second_clock::local_time();
            }
            k++;
        }


        if(settings->output.save_properties_evolution)
        {
            Log() << "Saving Properties Evolution\n";
            database.StorePropertiesEvolution(*settings, *prop);
        }
        if(settings->output.save_final_configuration)
        {
            Log() << "Saving Final Lattice\n";
            database.StoreFinalLattice(*settings, *lattice);
        }
        if(settings->output.save_final_properties)
        {
            Log() << "Saving Final Properties\n";
            PRE79MeanProperties pp(prop.get(), H.get());
            database.StoreFinalProperties(*settings, pp);
        }

        Log() << "Done\n";
        Log() << "Mean EPM: " << prop->TemporalMeanEnergyPerMolecule().Print() << std::endl;
        Log() << "Specific Heat: " << prop->SpecificHeat().Print() << std::endl;
        return *lattice;

    }
    virtual loggable_proxy<ILoggable> & Log()
    {
        ILoggable::Log() << "Thread: " << omp_get_thread_num() << "/" << omp_get_num_threads() << ": ";
        return ILoggable::Log();
    }
    /*

    ///for* parallel tempering
    double SwapTemperature(const double & t) {
        double oldT = H->GetTemperature();
        double oldL = H->GetLambda();
        double oldTau = H->GetTau();
        double oldH = H->GetH();

        delete H;
        H = new PRE79StandardHamiltonian(t,oldL,oldTau,oldH);
        return oldT;
    }
    */

    const Lattice * GetLattice() const
    {
        return lattice.get();
    }
    const PRE79StandardProperties * GetProperties() const
    {
        return prop.get();
    }
    const PRE79StandardProperties * GetThermalizationProperties() const
    {
        return thermalprops.get();
    }
    int GetNProductions() const
    {
        return productions.size();
    }
    const PRE79Production * GetProduction(const int & i) const
    {
        return productions[i].get();
    }
    const LatticeSimulation * GetSimulation() const
    {
        return simulation.get();
    }
    const LatticeSimulation * GetThermalization() const
    {
        return thermalization.get();
    }

};


#endif  /* _PRE79SIMULATION_H */

